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'INTRODUCTION 
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The  objectives  of  the  research  program  are  to  assess  how  well 
current  numerical  methods  and  turbulence  models  can  combine  to 
predict  the  detailed  flow  and  heat  transfer  behaviour  in  flow 
around  a  180'  bend,  a  geometry  of  great  importance  in  heat 
exchanger  design  and  in  various  other  cooling-passage 
geometries.  To  allow  this  assessment  detailed  flow  and 
thermal  experimental  data  are  required  -  and  these  were  not 
available.  The  total  research  program  therefore  has  two 
roughly  equal  parts:  obtaining  the  experimental  data  and 

developing  and  applying  the  3vd imens ional  numerical 
calculation  scheme. 


4b  e“- 

The  research  has  been  conceived  and  executed  as  a  joint  effort 
with  Professor  J.A.C.  Humphrey  and  his  group  at  UC  Berkeley. 
Mutual  help  has  been  provided  by  the  two  groups  over  the 
spectrum  of  tasks  but  it  would  be  correct  to  say,  in 
measurements,  that  the  main  emphasis  at  Berkeley  has  been  on 
flow  dynamics  while  that  at  UMIST  has  been  on  convective  heat 
transfer  aspects.  On  the  computational  side,  the  code 
development  for  the  square  sectioned  duct  has  mainly  proceeded 
at  Berkeley  while  that  for  the  circular  sectioned  tube  is 
proceeding  at  UMIST  '  (though  the  early  development  had  been 
undertaken  by  Professor  Humphrey  in  the  course  of  his  doctoral 
research) . 

Th£  present  report  describes  work  completed  at  UMIST  in  the 
past  year  in  furthering  the  accomplishment  of  the  research 
objectives.  Section  2  summarizes  the  computational  work  while 
progress  in  the  experiments  in  presented  in  Section  3. 
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2.  COMPUTATIONAL  RESEARCH  ON  FLOW  AROUND  180°  BENDS 

2.1  Tubes  of  Circular  Cross-Section 


2.1.1  The  starting  point  and  the  path 

The  starting  point  for  this  work  has  been  the 

three-dimensional  elliptic  flow  code  TOROID.  This  program  was 
developed  by  Professor  J.A.C.  Humphrey  while  a  Ph.D.  student 
at  Imperial  College,  London  flj  .  It  is  a  finite-volume 
discretization  of  the  steady-flow  Navier-Stokes  equations 
expressed  in  toroidal  coordinates.  It  follows  the  structure 
and  strategy  of  the  TEACH  family  of  computer  programs 
developed  at  Imperial  College  during  the  1970's  through  the 
guidance  and  overall  coordination  of  Dr.  A.D.  Gosman. 

Our  aims  in  relation  to  TOROID  are  twofold.  The  first  is  to 
improve  the  numerics  so  that  the  code  can  generate,  at 
moderate  computing  cost,  what  are  sensibly  accurate  solutions 
of  the  equations  of  motion  for  the  flow  around  a  180  bend. 
The  second  is  to  embed  a  high-level  turbulence  model  (a 
so-called  "Algebraic  Stress  Model"  -  ASM)  to  allow  reliable 
predictions  of  flow  pattern  and  local  heat  transfer  rates  over 
the  turbulent  flow  regime  for  arbitrary  flow  entry  conditions. 
The  work  completed  so  far  has  been  directed  at  the  first 
objective.  The  adaptations  we  wished  to  introduce  were: 

a.  The  incorporation  of  quadratic  upstream  differencing  of 
convective  transport  (QUICK)  in  place  of  hybrid 
upwind  differencing. 

b.  The  improvement  of  boundary-condition  treatments  near  the 
axis  in  order  to  remove  the  need  for  a  fine  grid  in  the 
centre  of  the  flow. 

c.  The  conversion  of  the  code  from  elliptic  to  "semi- 
elliptic"  in  order  to  allow  the  use  of  fairly  fine  meshes. 

d.  The  inclusion  of  the  SIMPLER  algorithm  for  handling  the 
pressure-continuity  connection  in  place  of  the  SIMPLE 
scheme  of  the  original  [2]. 

Items  a.  -  c.  have  been  completed  and  task  d.  is  well 
advanced;  technical  details  of  the  work  are  given  in  sections 

2.1.2  -  2.1.4. 
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The  further  steps  in  the  code  development  will  then  be 

e.  Inclusion  of  ka,c  Boussinesq  viscosity  model  of  turbulence 
(BVM) . 

f.  Inclusion  of  wall  functions  developed  earlier  in  the 
project . 

g.  Incorporation  of  a  uniform-property  ASM. 

h.  Inclusion  of  quasi-buoyant  effects  (associated  with 
centrifugal  pressure  gradients)  as  a  result  of  which 

the  heavier  eddies  will  be  flung  preferentially  to  the 
outside  of  the  bend. 

Testing  and  comparison  with  experiment  will  of  course 
accompany  each  stage  of  development. 

2.1.2  The  Describing  Equations 

Figure  2.1  defines  the  toroidal  coordinate  system  used  to 
describe  the  flow  field.  In  these  coordinates  the  continuity, 
momentum  and  energy  equations  for  a  steady,  uniform-density, 
high  Reynolds,  number  viscous  flow  may,  following  Aris  [3]  , 

for  example,  be  written: 

Continuity 

I*  rcU  +  fe  rrcV  +  Te  rW  "  0  {2a) 

where  rc  is  the  local  radius  •  R  *  rcos$  ,  and  U,  V  and  W 
denote  velocity  components  in  the  ♦,  r  and  0  directions 
respectively. 


Momentum  and  Energy  Equations 


p(C(*j  ♦  S  (*))  -  D(*)  ♦  Sn(*)  ♦  S  (♦)  (2.2) 

c  up 

where  $  stands  for  any  of  the  velocity  components  or  the 
temperature  T  and  the  operators  C(4>)  and  D(<0  stand  for: 

CM  -  -ji-  t  (rcU*)  +  ^  <rrcV*>  ♦  ^  (rW*)l  (2.3) 

c 

«♦>-  rrTr  ^  r  I?  (rcu  !♦>*!?  *rrcu  Tr  ^ 
c  if 


(2.4) 
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where  Pr,),  is  unity  for  the  velocity  components  and  denotes 
the  Prandtl  number  in  the  energy  equation.  The  functions 
Sc  (*)  ,  Sp('I-)  and  SD  (<i,)  are  given  in  Table  2.1.  It  will  be 
noted  that  viscous  terms  associated  with  O-direction 
second  derivatives  have  been  dropped,  a  simplification  that  is 
allowable  in  high  Reynolds  number  flows. 

2.1.3  Discretization 

Following  the  strategy  of  the  TEACH  codes  and  other 
finite-volume  primitive-variable  procedures,  the  differential 
equations  presented  in  Section  2.1.2  are  discretized  by 
integrating  them  over  small  contiguous  control  volumes  which 
collectively  cover  the  solution  domain.  To  each  control 
volume  is  attached  a  discrete  value  of  the  dependent  variable. 
A  staggered  arrangement  of  grid  nodes  is  adopted  as  indicated 
in  figure  2.2.  This  choice  minimizes  the  amount  of 
interpolation  that  has  to  be  done  and  improves  stability. 

In  evaluating  diffusion  processes  <D(<P))  a  linear  variation  of 
the  dependent  variables  between  adjacent  nodes  is  assumed 
while  the  source  terms  ( S  ( <^-)  ,  ^  (<!» )  and  Sc  (*))  are  mainly 
treated  as  uniform  over  any  particular  control  volume.  This 
is  the  standard  TEACH  practice. 

In  the  original  TOROID  code,  as  provided,  convective  transport 
was  handled  using  either  an  upwind  or  a  central-difference 
interpolation  depending  on  whether  the  modulus  of  the  cell 
Peclet  number  was  greater  than  or  less  than  2.  The  upwind 
approximation,  though  admirably  stable,  is  well  known  to  give 
rise  to  severe  false  diffusion  if  flow  lines  cut  at 
appreciable  angles  to  the  cell  boundaries  (the  error  being 
proportional  to  the  sine  of  twice  the  angle  of  incidence) .  A 
preliminary  study  during  the  first  year  of  ONR  funding  (Han, 
Humphrey  and  Launder  [4])  showed  decisively  superior  results 
from  using  Leonard's  (5]  quadratic  upwind  approximation 
(QUICK) .  The  principal  test  case  in  that  study  was  the  flow 
in  a  driven  cavity  which  produced  a  rather  similar  flow 
pattern  to  that  which  was  expected  in  the  cross-sectional 
plane  of  the  curved  tube.  At  a  cavity  Reynolds  number  of  1000 
roughly  the  same  numerical  accuracy  was  achieved  with  QUICK 
using  one  quarter  as  many  nodes  as  with  the  upwind/central 
approximation.  Accordingly,  at  the  outset  QUICK  was 
incorporated  into  the  TOROID  code.  In  devising  an 
approximation  for  the  value  of  at  a  point  e  midway  between 
nodes  i  and  i  +  1  one  Interpolates  on  a  parabola  passing 
through  these  two  points  and  through  the  next  point  on  the 
upwind  side  of  e.  Thus 

forUe>0  -  1<VW  *  J  <*l-l  *  ♦iel  -  2V 
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Table  2.1 


Source  terms  for  dependent  variables 


Toroidal  co-ordinate  system 


I  Uniform-9  surface 

-jj-  Control  volume  boundary 
for  P 

-II-  Control  volume  for  V 
••||*  Control  volume  for  U 
-  —  Control  volume  for  W 
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Uniform- 0  surface 


Fig. 


2,2 


Staggered  nodes  and  associated  control  volumes 
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Ue  <  0  *e  -  ««-£  ♦  *i+1>  -  i  <*i  +  X-i+2  -  2+i> 


In  fact,  for  the  purpose  of  improving  stability  the  formulae 
are  reorganized  as  follows 

ue  >  0  *e  “  i  (5*i  *  ^i+l*  "  3(Vl  "  V 

(2.5) 

Ue  <  °  *e  -  i  (3*.  *  5*.+1)  -  |(*.+2  -  *-+1> 

with  the  terms  in  the  right-hand  parenthesis  being  treated  as 
'source'  terms,  and  evaluated  using  previous-iteration  values. 

The  above  treatment  is  applied  in  TOROID  only  for  the  U  and  V 
component.  The  marching  character  of  the  solution  in  the 
general  flow  direction  means  that  an  upwind  approximation  must 
be  retained  for-  the  W  component.  In  any  event,  an  upwind 
approximation  is  much  more  satisfactory  for  the  streamwise 
than  the  cross-stream  components  for  the  flow  crossing  the 
upstream  and  downstream  faces  does  so  nearly  orthogonally  to 
the  control  volume  face  over  most  of  the  cross  section. 
Figure  2.3  compares  streamwise  velocity  distributions  obtained 
at  35°  from  the  ent^  to  a  circular  bend  for  a  simulation 
covering  the  first  70°^  of  arc.  The  entering  velocity  profile 
was  the  parabolic  distribution  associated  with  fully  developed 
flow  in  a  straight  tube.  Two  different  node  densities  are 
employed  (the  indicated  numbers  of  nodes  refer  to  the  8  ,  ♦ 
and  r  directions  respectively) .  There  is  a  moderate 
difference  between  the  profiles  for  QUICK  and  hybrid  (i.e. 
upwind/central)  treatments  with  the  finer  meshes.  Using  the 
coarser  mesh,  QUICK  leads  to  a  distribution  of  velocity 
similar  to  the  fine-grid  hybrid  scheme.  Having  regard  for  the 
earlier  calculations  Of  ref  (4)  there  was  little  doubt  that 
the  QUICK  results  were  the  more  accurate.  Even  for  that 

scheme,  however,  one  really  needs  a  finer  mesh  than  any  used 
in  the  above  tests  to  reduce  numerical  errors  to  insignificant 
levels.  A  modification  that  allows  the  desired  mesh 
refinement  is  summarised  in  section  2.1.5. 
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Flg.  2.3 


Coarse  gu.d  computations  comparing  HYBRID 
and  QUICK  differencing 
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2.1.4  Boundary  Conditions 


We  consider  flow  in  a  round  sectioned  tube  bent  in  a  circular 
arc.  The  flow  pattern  is  assumed  symmetric  about  the  plane 
<ji  =  0,  «  passing  through  the  centre  of  the  tube.  The  gradient 
of  V,  W,  p  and  T  normal  to  the  plane  is  therefore  zero  as  is 
likewise  the  value  of  U  (the  velocity  normal  to  the  plane)  . 
Around  the  tube  surface  all  velocity  components  are  set  to 

zero . 

There  has  been  some  controversy  in  the  literature  over  the 
handling  of  boundary  conditions  at  the  pipe  centre  (r  =  0). 
Provided  the  grid  is  made  fine  enough  the  assumptions  made 
about  the  boundary  values  become  unimportant  (in  illustration, 
one  can  imagine  that  a  fine  but  rigid  wire  bent  to  follow 

exactly  the  locus  of  the  pipe  centre  would  hardly  alter  the 
flow  pattern  from  that  found  with  the  wire  absent).  However, 
the  central  region  is  not  one  where  one  would  wish  to 

concentrate  a  fine  mesh  and  the  coarser  the  grid  the  more  care 
needs  to  be  taken. 

In  the  original  version  of  TOROID  the  radial  gradient  of  all 
dependent  variables  was  set  to  zero  at  the  axis.  With  only 
nine  radial  nodes  this  led  to  unphysical  features  in  the 

solution  and  accordingly  the  treatment  was  reformulated. 
Figure  2.4  shows  the  locations  of  nodes  in  the  r^v,  $  plane 
near  r  =  0.  A  central  problem  is  that  at  the  central  point 
many  grid  nodes  are  present  (each  corresponding  to  a  different 
value  of  $) .  For  the  streamwise  velocity  component  W  the  same 
value  was  assigned  to  all  the  central  nodes  being  the  average 
value  of  W  over  W(I,2): 

NI-1 

W(I,1)  =  t  W( I , 2) / (NI-2) 

2 

The  U  and  V  components  cannot,  however,  take  the  same  value  at 
the  centre.  Instead  we  required  that  the  resultant  of  each  U 
and  V  velocity  component  must  all  produce  the  same  velocity 
vector.  Now  the  resultant  velocity  must  lie  along  the 
symmetry  axis:  its  value  is  taken  as  the  average  of  the 

radial  velocity  components  on  the  symmetry  axis  either  side  of 
the  centre  node 

Vres.=0.5  (V(2, 3)  -  V(NI-1,3)) 

The  U  and  V  velocity  components  at  the  centre  are  then 
obtained  as 


I 
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Fig.  2.4 


Boundary  conditions  at  tube  centre 
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U(I,1)  =  V  sin(«u(I)) 
res  u 

V(  I  >  1 )  =  Vres  COS(HI)) 

Finally,  the  ring  of  nodes  V(I,2)  is  assigned  by  linearly 
interpolating  between  V(I,1)  and  V(I,3).  (This  last  step  is 
not  an  essential  element  of  our  scheme:  one  could 

alternatively  obtain  V(I,2)  by  solving  the  radial  momentum 
equation  over  the  wedge  shaped  region  shown  in  figure  2.4  as 
is  currently  being  done  by  one  of  Professor  Humphrey's 
students  at  UC  Berkeley) . 

Figure  2.5  compares  profiles  of  radial  velocity  across  the 
plane  of  symmetry  for  the  original  and  new  treatments.  Again 
the  computations  wereQ  those  of  a  70°  bend  and  the  profiles 
shown  are  those  at  45°.  The  new  version  has  entirely  removed 
the  unphysical  discontinuities  in  velocities  that  arose  with 
the  original. 


2.1.5  The  Semi-Elliptic  Procedure 

A  semi-elliptic  solving  procedure  is  one  where  in  one  of  the 
coordinate  directions,  for  all  dependent  variables  except  the 
pressure;  information  is  held  to  be  transmitted  only  one  way, 
i.e.  from  "upstream"  to  "downstream".  The  practical  outcome 
of  this  supposition  is  that  these  variables  do  not  require 
storing  over  the  whole  solution  domain;  instead  just  two 
planes  of  storage  are  provided  for  each  variable 
(corresponding  to  the  current  plane  and  that  immediately 
upstream  thereof)  and  the  computation  "marches"  from  upstream 
to  downstream  successively  overwriting  upstream  values  held  in 
an  array  by  downstream  values  on  the  completion  of  each  new 
forward  step.  The  pressure  must  be  treated  differently  from 
other  variables  for  (in  subsonic  flow)  it  is  intrinsically 
affected  by  the  flow  in  all  directions.  Indeed,  for  the  flow 
in  a  pipe  bend,  it  is  the  upstream  transmission  of  pressure 
information  from  downstream  that  "warns"  the  flow  that  it  is 
approaching  the  bend  and  will  have  to  change  direction 
accordingly.  The  pressure,  therefore,  must  be  held  over  the 
entire  domain  and  inevitably  the  solution  -  not  to  just  the 
pressure  but,  by  virtue  of  the  coupling  of  the  variables,  to 
the  velocity  components  too  -  must  be  obtained  by  making  many 
sweeps  over  the  domain  successively  refining  the  pressure 
field  at  each  pass. 

Thus,  a  semi-elliptic  procedure  does  not  offer  a  faster  route 
to  solution.  By  drastically  reducing  the  storage 
requirements,  however,  it  allows,  for  a  given  computer 
installation,  a  far 


OLD  TREATMENT 


NEW  TREATMENT 


2.5 


Comparison  of  axial  and  radial  velocity 
profiles  with  new  and  centre  boundary 
treatments 
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greater  density  of  grid  nodes  than  with  a  fully  elliptic 
scheme  and  thereby  brings  within  the  sights  of  computational 
feasibility  a  vast  range  of  complex  engineering  flows.  The 
180°  tube  bend  is  one  flow  which,  with  a  CDC  7600  computer, 
becomes  definitely  computable  with  a  semi -el 1 i ptic  scheme. 

At  high  Reynolds  numbers  the  only  restriction  on  adopting  a 
semi-elliptic  treatment  is  that  there  should  be  no  reverse 
flow  regions  in  the  direction  of  marching.  This  prevents 
consideration  of  flow  around  very  tight  bends  where  separation 
occurs. 

The  detailed  strategy  of  a  semi-elliptic  solving  procedure  has 
been  set  out  by  Pratap  and  Spalding  [6]  and  others.  A 
somewhat  different  sequencing  is  required  than  with  a  fully 
elliptic  solver  due  to  the  fact  that  no  velocity  field 
information  is  available  from  the  previous  streamwise  sweep 
over  the  flow.  Figure  2.6  indicates  the  steps.  The  diagrams 
show  N,  W  and  P  nodes  lying  on  a  constant- $  surface. 
(Although  the  cells  are  shown  as  rectilinear  a  small  curvature 
-  about  2°  of  arc  ~  is  present  in  the  6  direction)  .  The  nodes 
are  shown  by  arrows  (for  velocity)  or  circles  (pressure); 
solid  symbols  indicate  that  the  nodal  value  has  been 
calculated  while  for  open  symbols  no  value  is  .available.  The 
half-filled  pressure  nodes  contain  the  value  of  pressure  from 
the  previous  sweep  over  the  flow  domain.  At  a  given  step  the 
sequence  begins  by  computing  current  plane  values  of  V  and  u 
(not  shown)  by  solving  the  respective  momentum  equations  using 
the  pressure  from  the  previous  sweep  and  values  for  velocity 
from  the  upstream  station.  The  streamwise  momentum  equation 
is  next  solved  to  give  the  W  velocity  (displaced  a  half  cell 
downstream)  using  "new"  values  of  V  and  U  in  the  convection 
terms  and  presuming  a  cell  mass  balance  in  obtaining  the  rate 
of  momentum  outflow  at  the  downstream  face  of  the  control 
volume . 

The  final  step  is  the  adjustment  of  the  pressure  field  by  the 
application  of  the  continuity  equation  to  the  current-plane 
pressure  cells.  Perturbations  in  the  pressure  and  U  and  V 
velocities  at  the  current  plane  and  perturbations  in  pressure 
(only)  at  the  upstream  plane  are  introduced  to  reduce  mass- 
balance  errors.  Finally,  a  bulk  correction  to  the  downstream 
pressure  is  applied  to  achieve  global  mass  conservation  over 
the  tube  cross  section. 

The  same  sequence  of  operations  is  then  applied  at  the  next 
step,  and  the  process  successively  repeated  until  the  whole 
domain  has  been  covered.  The  calculation  then  begins  a 
further  sweep,  starting  at  the  upstream  plane  and  proceeding 
step-by-step  over  the  flow  domain  using  the  updated  pressure 
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of  solution  with  semi-elliptic 


field  in  solving  the  momentum  equations  to  compute  the 
velocities. 

Because  the  computation  as  described  above  has  had  to  make 
extensive  use  of  "upstream*  values  instead  of  current-plane 
values  in  evaluating  sources  and  convection  coefficients,  even 
when  the  mass  residuals  have  been  reduced  to  negligibly  small 
levels  the  pressure  and  velocity  fields  will  still  differ 
somewhat  from  those  obtained  by  a  fully  elliptic  solution.  It 
is  therefore  necessary,  as  the  calculation  approaches  its 
(apparent)  converged  solution  to  introduce  iteration  on  the 
velocity  components  at  each  step.  That  fs  to  say,  when 
"current-plane"  values  have  been  obtained  the  velocity 
equations  are  re-solved  using  current  values  as  appropriate  in 
re-forming  the  coefficients  and  source  terms*.  Figure-  2.7 
shows  the  very  substantial  change  in  the  secondary  velocity 
that  results  from  this  practice.  Due  to  the  non-linearity  of 
the  equations  three  iterations  are  needed  to  achieve  sensibly 
asymptotic  values  (though  for  industrial  calculations  or  for 
turbulent  flow  the  additional  computing  cost  would  probably 
not  justify  going  beyond  the  second  iteration) . 


2.1.6  Application  to  Flow  around  a  180°  Bend 

Agrawal  et  al  [7]  have  provided  a  laser  Doppler  study  of 
laminar  flows  around  a  180°  bend  over  a  range  of  Reynolds 
numbers  and  Dean  numbers.  From  these  the  test  at  a  Dean 
number  of  183  has  been  selected  for  computer  simulation 
with  the  TOROID  program.  In  the  computations  a  uniform 
streamwise  velocity  and  static  pressure  have  been  assumed. 
The  former  corresponded  closely  with  the  inlet  velocity 
profile  reported  by  the  experimenters  though  inevitably  there 
was  a  thin  boundary  layer  present  whose  possible  effects 
cannot  be  entirely  discounted. 


♦Pratap  and  Spalding  [6]  apparently  omitted  this  corrective 
step.  That  they  should  have  still  obtained  satisfactory 
agreement  with  experiment  is  presumably  because  the  pipe  bend 
considered  in  their  study  was  much  less  severe  than  that 
which  we  are  tackling. 


upstream  coefficients  ' 

I 

Fi9*  2.7  Comparison  of  semi-elliptic  field  with  elliptic 

solution 
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A  uniform  20x20  mesh  has  been  used  to  map  the  semicircular 
cross  sectional  plane  with  100  equispaced  streamwise  sections. 
Successive  sweeps  were  made  until  the  mass  and  momentum 
residual  errors  summed  over  all  the  sections  were  reduced 
below  10~3  0f  the  entering  mass  and  momentum  flows.  At  exit 
zero-gradient  boundary  conditions  were  applied.  This  weak 
constraint,  while  not  particularly  accurate,  would  not  be 
expected  to  contaminate  the  solution  more  than  a  few  tens  of 
degrees  inboard  from  the  exit  plane.  As  it  happened  the  last 
station  at  which  Agrawal  report  velocity  profiles  was  at 
145  degrees,  at  which  position  it  seems  unlikely  there  would 
be  significant  effect  on  the  numerical  solution  due  to 
uncertainties  in  the  outlet  pressure  field. 

Figure  2.8  compares  the  measured  and  computed  streamwise 
velocity  profiles  along  five  lines  across  the  pipe  cross 
section  and  at  five  streamwise  positions.  The  figures  of  the 
measured  distributions  are  taken  directly  from  the  published 
paper;  the  computer  plotted  numerical  solutions  are  shown  in 
mirror  image  and  to  a  slightly  different  scale.  It  is 

nevertheless  clear  that  there  is  a  close  degree  of 
correspondence  between  the  measured  and  computed  behaviour. 
Perhaps  the  most  significant  difference  between  the  two  is 
seen  at  the  first  station.  The  computed  velocity  along  the 
line  running  closest  to  the  wall  is  markedly  higher  than 
measured  while  along  the  other  lines  agreement  is  very  good. 
The  discrepancy  almost  certainly  arises  from  neglecting  a  thin 
boundary  layer  at  entry  to  the  bend.  Since  it  is  the  velocity 
defect  of  this  boundary  layer  that  provides  the  source  term 
driving  the  secondary  flow  around  the  outer  periphery  it  is 
consistent  that  we  notice  at  1  >  160°  (L/R  *  19.53) 
the  computed  velocity  profiles  do  not  exhibit  quite  as  much 
distortion  as  the  measurements,  an  indication  that  the 
calculated  secondary  flow  is  slightly  less  than  in  the 
experiment.  Unfortunately,  the  experimenters  do  not  report 
secondary  velocities  at  the  same  Dean  number  as  for  the  axial 
flow.  Comparisons  of  the  present  calculations  (for  a.  Dean 
numberQof  183)  with  experiments  at  a  Dean  number  of  138  at 
6  -  30  shown  in  Figure  2.9  do  display  broadly  similar 
distributions  though,  as  would  be  expected,  the  magnitude  of 
the  computed  flow  is  somewhat  greater  than  measured. 


STREAMWISE  VELOCITY  PROFILES 
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2.2  Rectangular  Sectioned  Duct 


2.2.1  Summary  of  computational  work 


The  program  at  UMIST  has  been  the  beneficiary  of  the 
semi-eil iptic  code  developed  at  Berkeley.  The  recent  report 
by  Humphrey  et  al.  [1]  describes  this  solving  procedure  and 
its  application  to  the  dynamics  of  flow  around  a  180°  bend  for 
a  duct  of  square  cross-section. 


The  UMIST  effort  has  consisted  of:  getting  the  UCB  code 
mounted  on  the  local  computer  system  (unfortunately  this 
proved  a  major  task  due  to  tape-read  problems  )  ;  incorporating 
new  wall  resistance  laws;  adding  the  energy  equation  to  allow 
the  thermal  field  to  be  calculated;  including  diagnostics  to 
give  a  more  complete  picture  of  residual  errors  in  solving  the 
'difference'  equations;  the  introduction  of  "iterations"  on 
variables  other  than  pressure  (as  discussed  for  the  toroidal 
duct  in  section  2.1);  developing  plotting  codes  for  printing 
the  results;  and,  finally,  making  the  computations  themselves. 

The  changes  made  to  the  code  are  discussed  in  section  2.2.2 
while  section  2.2.3  presents  the  computational  results  and 
draws  comparison  with  those  obtained  at  Berkeley. 


2.2.2  Principal  changes  to  the  Berkeley  curved  duct  code 


a)  multiple  iterations 


Just  as  with  TOROID,  It  was  found  that  accurate  solution  of 
the  momentum  equations  required  iteration  of  the  dependent 
variables  at  each  plane.  This  practice  allows  source  terms 
and  convective  fluxes  which,  in  the  original  version,  had  of 
necessity  to  be  approximated  by  corresponding  values  at  the 
adjacent  upstream  plane,  to  be  re-evaluated  locally. 
Streamwise  variations  are  so  rapid  in  certain  regions  of  the 
flow  around  the  bend  that  a  considerably  finer  mesh  than  the 
3°  intervals  used  in  the  present  work  would  have  been  required 
to  make  the  use  of  upstream  values  fully  satisfactory. 


After  explorations  with  various  numbers  of  local  momentum 
iterations,  the  practice  of  3  iterations  per  plane  was  applied 
as  standard.  This  figure  could  however  have  been  dropped  to  2 
without  altering  the  computed  values  of  velocity  by  more  than 
1%  with  a  saving  of  'some  22%  in  computer  time.  (In  this 
context,  '1  iteration'  would  mean  the  original  practice  of 
using  upstream  coefficients  and  proceeding  wi thout  iteration.) 


I  lliven  the  weaknesses  in  the  present  physical  model  there  is  no 
I  Practical  benefit  in  solving  the  numerical  equations  to  the 
fndicated  accuracy.  In  this  research,  however,  it  is 
/iesirable  to  identify  the  origin  of  computational  error, 
i  whether  numerical  or  physical. 

1 15 )  Aids  to  Convergence 

The  multiple  iteration  practice  noted  above  proved  to  be 
highly  destabilizing  if  introduced  early  in  the  iteration 
sequence.  Accordingly,  only  when  the  pressure  field  was 
approaching  its  'converged*  solution  (without  iteration)  was 
iteration  on  'the  momentum  equations  introduced;  first  one 
iteration,  later  two  etc.  In  addition  formal 

upstream-downstream  under-relaxation,  as  introduced  in 
Pratap's  initial  semi-elliptic  computations, 
was  used  for  stability  in  the  initial  stages.  This  practice 
j  (as  noted  earlier  with  respect  to  TOROID)  had  to  be  switched 

out  before  the  end  of  the  convergence  to  avoid  errors. 

As  with  TOROID,  a  diagnosis  of  the  accuracy  of  convergence  of 
solution  of  the  different  equations  over  the  whole  domain  was 
added  for  all  dependent  variables.  With  a  'zero'  initial  field 
for  pressure,  45  sweeps  around  the  bend  were  needed  to  reduce 
the  global  normalised  residual  errors  below  1%  for  a  12x22 
grid  in  the  cross-stream  plane  and  103  streamwise  planes  (of 
which  60  were  in  the  180°  bend  section)  .  The  central 
processor  time  required  for  such  a  calculation  for  a  CDC  7600 
computer  was  1280  seconds  (decimal) . 

I  2.2.3  The  computed  behaviour  for  the  flow  in  a  square 

)  sectioned  duct  around  a  180°  bend. 

The  present  section  discusses  computational  results  for  the 
;  test  geometry  studied  in  the  ONR-funded  work,  a  ratio  of  mean 

;  bend  radius  to  hydraulic  diameter  of  3.35:1.  The  Reynolds 

i  number,  56000,  is  that  used  by  Professor  Humphrey's  team  in 

their  measurements  of  the  same  flow.  The  present  results  have 
,  all  been  obtained  with  the  k- e  Boussinesq  viscosity  model 

|  j  using  standard  coefficients.  The  interest  at  this  level  of 

|  I  modelling  is:  to  observe  the  sensitivity  of  the  results  to 

|  different  wall-boundary  treatments,  to  compare  our  numerical 

|  I  results  with  those  already  obtained  by  Professor  Humphrey's 

|  group  for  the  same  physical  model,  and  to  observe  the 
predicted  distribution  of  heat  transfer  coefficients.  These 
,  aspects  are  considered  in  turn. 
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An  impression  of  the  effect  of  the  new  wall  treatment  on  the 
veolcity  field  emerges  from  figure  2.10a.  The  135°  position 
was  chosen  as  this  displayed  the  maximum  difference  between 
the  two  treatments,  i.e.  the  standard  (I.C.)  treatment  and  the 
modified  (R.J.)  version.  For  the  present  flow  the  most 
significant  difference  in  the  two  approaches  is  that,  in  the 
•RJ'  scheme,  the  wall  friction  opposing  the  secondary  motion 
is  obtained  independently  of  the  streamwise  velocity  by 
performing  the  integration 

v  =/{  Ty*  > 

J  Weff 

between  the  wall  and  the  first  node.  In  contrast,  the  IC 
wall  treatment  assumes  that  the  resultant  near-wall  velocity 
parallel  to  the  surface  obeys  the  usual  logarithmic  law  (the 
streamwise  and  cross-stream  components  are  then  obtained  by 
resolving  appropriately).  There  is  very  little  difference 
between  the  velocity  fields  generated  by  the  two  approaches. 
The  effect  on  the  resultant  stress  and  heat  flux  fields  is 
greater  but  still  relatively  small  (up  to  10%).  The  result 
does  not  indicate'  that  the  flow  is  essentially  independent  of 
the  nature  of  the  wall  law  for  the  secondary  motion.  When  the 
wall  stress  opposing  the  secondary  flow  is  set  to  zero  it  is 
seen  from  figure  2.10b  that  quite  substantial  modifications  of 
the  streamwise  and  axial  velocity  fields  occur.  These  results 
indicate  that  we  may  have  some  difficulty  in  sifting  out 
whether  errors  due  to  the  physical  model  arise  from  the  wall 
treatment  or  the  interior  turbulence  model.  The  UMIST  heat 
transfer  data  will  be  very  helpful  in  this  regard,  because, 
they  will  be  much  more  sensitive  to  the  wall  treatment  than  is 
the  interior  velocity. 

Comparisons  between  the  present  computations  of  the  velocity 
field  and  those  obtained  by  Professor  Humphrey's  team  are 
shown  in  figure  2.11  at  45  and  90  .  Results  obtained  from 
using  both  QUICK  and  HYBRID  treatments  of  convection  are 
presented  though  the  former  has  been  established  as  the  more 
accurate  scheme  (4  ,  8].  There  is  broadly  the  same  behaviour 
shown  by  the  independent  computations  made  at  the  two 
institutions:  our  results  confirm  the  prior  discovery  by  our 
Berkeley  colleagues  that,  with  the  k-e  Boussinesq  model,  both 
the  streamwise  and  secondary,  flow  patterns  are  in  substantial 
disagreement  with  experiment  from  45°  onwards. 


O 
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There  are,  nevertheless  some  differences  between  the  computed 
results  obtained  at  the  two  institutions  that  are  more  serious 
in  the  case  of  the  numerically  more  accurate  QUICK  scheme(for 
clarity  the  Berkeley  results  with  HYBRID  are  omitted).  There 
are  several  factors  possibly  contributing  to  the  differences. 
To  mention  but  two,  the  UMIST  results  used  a  coarser  grid  than 
the  Berkeley  computations  (to  allow  us  to  exploit  the  local 
computer  centre  limits  on  job  sizes);  on  the  other  hand  the 
present  results  used  current-plane,  iterated  values  of  all 
dependent  variables  whereas  .  it  appears  that  the  Berkeley 
results  may  (in  common  with  earlier  semi-elliptic  treatments) 
have  used  upstream  values.  Clearly  the  origin  of  the 
differences  needs  identifying  and  an  agreed  definitive  set  of 
results  produced.  It  is  our  view  that  the  discovery  of 
variations  between  our  results  provides  support  for  there  to 
be  planned  overlap  between  research  projects  a  feature  that, 
though  incorporated  from  the  beginning  in  this  Berkeley  - 
UMIST  collaboration,  is  rarely  built  into  research  programs. 

The  computed  heat  transfer  pattern  in  the  flow  around  the  bend 
is  shown  in  figure  2.12.  Experimental  data  are  not  yet 
available  with  which  to  draw  direct  comparison  though,  given 
the  rather  poor  quality  of  agreement  of  the  velocity  field,  it 
is  certain  that  there  will  be  some  significant  differences. 
Nevertheless  the  broad  trends  in  the  observed  heat  transfer 
rates  are  probably  not  in  question.  The  level  of  the  Nusselt 
number  is  strongly  affected  by  the  secondary  flow.  The 
migration  of  fluid  radially  outward  near  the  mid-plane 
produces  slow  moving  relatively  sluggish  fluid  near  the  inner 
wall  and  associated  low-levels  of  heat  transfer  coefficient. 
Correspondingly,  fluid  impingement  on  the  outer  wall  causes 
steep  gradients  of  streamwise  velocity  and  high  Nusselt 
numbers.  The  ratio  of  maximum  to  minimum  values  of  Nu  exceeds 
2:1  over  the  second  90°  of  the  bend  and  downstream  therefrom 
more  than  10  hydraulic  diameters  (the  extent  of  the  present 
computations) .  It  is  likely  that  the  actual  level  of  Nusselt 
number  would  display  rather  greater  variations  than  this 
between  the  inner  and  outer  surfaces  of  the  bend;  for,  the 
exact  generation  processes  associated  with  mean  strain  will 
tend  to  augment  turbulence  activity  near  the  concave  outer 
wall  while  causing  damping  near  the  convex  inner  surface.  The 
algebraic  stress  model  we  shall  be  incorporating  in  the  next 
phase  of  our  modelling  work  accounts  exactly  for  these 
generation  processes  but  the  simpler  Boussinesq  k-e  model 
used  here  does  not. 


2*  (X;  Y)  /D 


45.00  DEGREES 


2*  (X;  Y)  /Dm 


90.00  DEGREES 


0.2  0.4 

0.6  0.8  1.0  1.2  1.4  1.6 

1.8 

2*  (X,  Y)  /Dh 

Z=  135.00  DEGREES 

CONCAVE 

RE=  56690.0' 

CONVEX 

'  QUICK  RJ  VT 

BOTTOM  (FULL  LENGTH) 

?ig.  2.12 

Computed  distribution  of  Nusselt 

number 

in  flow  around  180°  bend  (square 
duct) 

fn>  rM  o n°  ns° 

section 

-14b- 


2*(X;Y)/DM  Z=  5.00  Dh  DOWNSTREAM 

n  n 


2x(X;Y)/Dh  Z=  10.23  Dh  DOWNSTREAM 

—  CONCAVE  RE=  56690.0 

—  CONVEX  QUICK  RJ  WT 

—  BOTTOM  (FULL  LENGTH) 


(d)  180°  (e)  5.0Dh  downstream  (£)  10Dh  downstream 


We  note  finally  that  the  Nusselt  numbers  along  the  plane 
sidewalls  are  fairly  uniform  in  the  bend  itself  and 
intermediate  in  magnitude  between  those  on  the  inner  and  outer 
curved  walls.  The  mean  Nusselt  number  given  by  the 
T ' ttus-Boelter  correlation  is  130  at  this  Reynolds  number 
whereas  the  mean  computed  Nusselt  number  over  the  second  90 
of  bend  is  in  the  range  180-195,  an  augmentation  of  some  40%. 
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1.  EXPERIMENTAL  PROGRAM 

.The  experimental  work  is  running  several  months  behind  the 
schedule  foreseen  in  the  last  annual  report.  One  of  the  two 
major  factors  contributing  to  the  delay  has  been  the  fact 
that  far  more  machining  work  to  flanges  has  proved  necessary 
than  in  the  original  design.  Although  the  technician  time 
involved  has  been  met  by  UMIST  (rather  than  the  contract 
budget)  there  has  inevitably  been  delay  because  no  additional 
manpower  was  available.  Moreover  the  emphasis  of  Richard 
Johnson's  effort  has  been  consciously  shifted  from  experiment 
to  computation,  in  order  that  we  could  obtain  predictions  of 
the  heat  transfer  distribution  prior  to  installing  the  wall 
thermocouples  in  the  bend.  After  having  obtained  the 
computations  of  Nusselt  number  described  in  the  previous 
section  it  did  indeed  appear  that  some  re-arrangement  from  the 
original  plan  was  desirable. 


The  following  summarizes  the  tasks  now  completed  and  those 
still  remaining  before  the  definitive  tests  can  be  undertaken. 

1.  The  straight  and  curved  sections  are  completely  built 
including  a  supporting  structure  for  the  apparatus. 

2.  The  fan  (in  suction  mode)  has  been  installed  connected  and 
is  fully  operational. 

Because  the  electrically  heated  film  would  quickly  burn 
out  if  the  fan  should  fail  for  any  reason  a  safety  cut-out 
to  the  heating  circuit  has  been  installed. 

3.  A  metering  pipe  with  orifice  plates,  has  been  added  at  the 
downstream  end  of  the  test  section  to  monitor  flow  rate. 
This  gives  a  more  accurate  measurement  with  better 
resolution  than  the  originally  planned  static  pressure 
tappings  in  the  contracting  section. 
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4.  Leak  tests  have  been  performed  showing  very  good  integrity 
with  small  leaks  easily  plugged. 

'5.  Electrical  voltage  regulating  and  monitoring  equipment 

has  been  built  and  is  ready  for  use  to  regulate  the  power 
applied  to  the  Intrex  heating  film,  including  regulation 
of  14  quarter-inch  wide  'tracks'  to  be  inscribed  in  the 
intrex  used  on  the  flat  surfaces  around  the  bend.  This 
special  treatment  is  necessary  to  ensure  a  uniform  heat 
flux  there  (a  discussion  of  tl^is  issue  is  presented  in  the 
1980/81  Annual  Report  [9]). 


6.  An  inlet  contraction  section  has  been  constructed  and 
installed  in  the  duct  to  provide  a  smooth,  uniform  inlet 
velocity  profile  without  separation  or  other 
unsteadiness. 

7.  Fifteen  thermocouples  at  each  of  8  axial  measuring 
.  stations  have  been  installed  into  the  wall  of  the 

Plexiglas  duct.  The  thermocouples,  will  be  in 
contact  with  the  back  (untreated)  side  of  the 
Intrex  with  silver-loaded  epoxy  positioned  at  their 
measuring  junctions  to  provide  accurate  temperature 
measurments. 

Tests  have  been  carried  out  to  measure  the  temperature 
drop  across  the  Intrex;  this  appears  to  be  in  the 
neighbourhood  of  0.1°C  or  less  for  the  heat  fluxes 
desired . 

8.  Flow  visualization  tests  have  been  performed 
using  atomized  paraffin  droplets  and  a  'thin  slice’ 
laser  beam  produced  by  placing  a  glass  tube  in  front 
of  the  beam.  The  'thin  slice'  of  light  produced 
was  approximately  3  mm  thick  and  wide  enough  to 
illuminate- the  entire  cross-section  of  the  duct. 

At  present  visualization  has  been  successful  only 
at  low  Reynolds  numbers  an  example  of  which  is 
shown  in  figure  3.1.  The  secondary  flow  is  clearly 
evident  (the  assymnetry  about  the  horizontal  plane 
apparent  in  this  figure  arose  from  blockage  from  the 
smoke  injection  since  at  the  time  smoke  had  been 
introduced  only  a  short  distance  upstream  of  the  bend)  . 
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9.  The  gold-coated  Intrex  film  is  about  to  be  installed 
A  good  deal  of  testing  of  the  Intrex  has  been  carried 
out  to  establish  its  resistivity  characteristics,  the 
techniques  for  mounting  it  on  the  Plexiglas,  the  best 
way  of  ensuring  good  electrical  contact  and  the  upper 
voltage  limit  permitted  across  the  sheet.  The  various 
shake  down  tests  have  confirmed  indications  that  the 
Intrex  is  very  suitable  for  our  testing  purposes  and  will 
be  easy  to  handle  and  install. 


The  remaining  tasks  before  the  testing  programme  can  begin 
are: 

(a)  finish  the  installation  and  checking  of  the  Intrex, 

(b)  final  leak  testing, 

(c)  checking  of  fan,  voltage-control  systems,  manometers. 


(d)  build  temperature  rake. 

The  initial  measurements  will  be  a  limited  mean  velocity 
survey  upstream  of  the  bend  following  which  the  wall 
conditions  in  the  first  lm  of  the  inlet  duct  will  be 
successively  adjusted  to  achieve  as  closely  as  possible 
conditions  identical  with  those  of  the  Berkeley  experiments. 
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